* implements a regulation as a productivity shock defined by four input values:
*
*   s_reg:          regulated sector set by defining s_reg("xxx") = YES;
*
*   ex_shock_share: share of the shock value associate with extant capital based
*                   production
*                   -1 = all capital is affected
*                   0 = only new capital is affected
*                   1 = only extant capital is affected
*
*   shock_type:     flag defining the input bias of the shock
*                   1 = hick neutral
*                   2 = nestor and pasurka vector based on air regulations
*                   3 = capital only
*                   4 = labor only
*                   5 = economy-wide shares of capital and labor
*
*   shock_size:     initial year value of the shock in billions
*
* the regulation is implemented as a per unit output cost increase for the
* affected sources as defined by s_reg and ex_shock_share. shock_size is used to
* define the per unit cost based on all benchmark year production from the
* regulated sector. this per unit cost is then applied to all affected
* production based on ex_shock_share. the shock's size therefore, scales with
* output across time and regions. the per unit cost is shared across inputs
* based on the method defined by shock_type.
*
* shock_size is interpreted such that if production in the regulated sector was
* leontief then the value of the additional inputs needed to produce the
* benchmark level of output (assuming all vintages are affected) would be equal
* to the shock_size. where in this case the inputs are valued at there social
* opportunity costs, that is market prices gross of taxes.
*

set
  s_reg(s)              regulated sectors;

parameter
  shock_type            flag defining the input bias of the shock
  shock_size            value of the shock in billions
  ex_shock_share        share of the shock associated with extant capital
  per_unit(s)           per unit cost from shock
  shock(t,r,s,*)        region - year - vintage specific shock size
  prod_share(*)         for non-hicks neutral shock types defines the input bias
  share(t,r,*,s,*)      share of the shock associated with each production input
  denom(t,r,s)          denominator in defining the shock share;

* there are no regulated sectors until one is set
s_reg(s) = NO;

* load the shock specification
$include examples/shock_definition.gms

* per unit output cost from the shock
per_unit(s_reg)$sum((tfirst,rr), y.l(tfirst,rr,s_reg)*y0(rr,s_reg)
                                +y_ex.l(tfirst,rr,s_reg)*y_ex0(rr,s_reg))
                       = shock_size/sum((tfirst,rr), y.l(tfirst,rr,s_reg)*y0(rr,s_reg)
                                                    +y_ex.l(tfirst,rr,s_reg)*y_ex0(rr,s_reg));

* the year, region, and vintage enigeering cost is determined under the
* assumption that shock_size defines a per unit output cost of the regulation
* for the affected sources
if (ex_shock_share eq -1,
  shock(t,r,s_reg,"new") = per_unit(s_reg)
                             *y.l(t,r,s_reg)*y0(r,s_reg)
                             *pr.l(t,r);
  shock(t,r,s_reg,"extant") = per_unit(s_reg)
                                *y_ex.l(t,r,s_reg)*y_ex0(r,s_reg)
                                *pr.l(t,r);
else
  shock(t,r,s_reg,"new") = per_unit(s_reg)
                             *(1-ex_shock_share)
                             *y.l(t,r,s_reg)*y0(r,s_reg)
                             *pr.l(t,r);
  shock(t,r,s_reg,"extant") = per_unit(s_reg)
                                *ex_shock_share
                                *y_ex.l(t,r,s_reg)*y_ex0(r,s_reg)
                                *pr.l(t,r);
);

* input shares based on the nestor and pasurka air regulation vector
if (shock_type eq 2,
  prod_share("ele") = 0.27;
  prod_share("con") = 0.06;
  prod_share("wpm") = 0.01;
  prod_share("chm") = 0.01;
  prod_share("ref") = 0.01;
  prod_share("prm") = 0.025;
  prod_share("cem") = 0.025;
  prod_share("cpu") = 0.006;
  prod_share("tem") = 0.001;
  prod_share("bom") = 0.003;
  prod_share("trn") = 0.01;
  prod_share("ttn") = 0.01;
  prod_share("srv") = 0.2;
  prod_share("l") = 0.16;
  prod_share("k") = 0.2;
);

* input shares based on a capital only shock
if (shock_type eq 3,
  prod_share("k") = 1;
);

* input shares based on a labor only shock
if (shock_type eq 4,
  prod_share("l") = 1;
);

* input shares based on economy-wide capital and labor shares
if (shock_type eq 5,
  prod_share("k") = sum((r,s), (1+tk0(r))*(kd0(r,s)+kd_ex0(r,s)))/
                    sum((r,s), (1+tk0(r))*(kd0(r,s)+kd_ex0(r,s))
                              +(1+tl0(r))*(ld0(r,s)+ld_ex0(r,s)));
  prod_share("l") = 1-prod_share("k");
);

* next we need to determine the share of the shock in a region region and year
* that is associated with each input into the production function for both new
* and extant capital - the approach to doing this is different depending on
* whether exogenous shares are used (capital only, labor only, or nestor and
* pasurka) or whether endogenous shares are being used, as in the hicks netural
* shock type
if (shock_type eq 1,

* input value over which hicks neutral shock is applied for new capital
  denom(t,r,s_reg) = (1+tk(t,r))*pr.l(t,r)*kd.l(t,r,s_reg)*kd0(r,s_reg)
                    +(1+tl(t,r))*pl.l(t,r)*ld.l(t,r,s_reg)*ld0(r,s_reg)
                    +sum(s,pa.l(t,r,s)*id.l(t,r,s,s_reg)*id0(r,s,s_reg));

* share of the shock for production with new capital
  share(t,r,"k",s_reg,"new")$denom(t,r,s_reg) = (1+tk(t,r))*pr.l(t,r)
                                                *kd.l(t,r,s_reg)*kd0(r,s_reg)
                                                /denom(t,r,s_reg);
  share(t,r,"l",s_reg,"new")$denom(t,r,s_reg) = (1+tl(t,r))*pl.l(t,r)
                                                *ld.l(t,r,s_reg)*ld0(r,s_reg)
                                                /denom(t,r,s_reg);
  share(t,r,s,s_reg,"new")$denom(t,r,s_reg) = pa.l(t,r,s)
                                              *id.l(t,r,s,s_reg)*id0(r,s,s_reg)
                                              /denom(t,r,s_reg);

* input value over which hicks neutral shock is applied for extant capital
  denom(t,r,s_reg) = (1+tk(t,r))*pr_ex.l(t,r,s_reg)*kd_ex.l(t,r,s_reg)*kd_ex0(r,s_reg)
                    +(1+tl(t,r))*pl.l(t,r)*ld_ex.l(t,r,s_reg)*ld_ex0(r,s_reg)
                    +sum(s,pa.l(t,r,s)*id_ex.l(t,r,s,s_reg)*id_ex0(r,s,s_reg));

* share of the shock for production with extant capital
  share(t,r,"k",s_reg,"extant")$denom(t,r,s_reg) = (1+tk(t,r))*pr_ex.l(t,r,s_reg)
                                                   *kd_ex.l(t,r,s_reg)*kd_ex0(r,s_reg)
                                                   /denom(t,r,s_reg);
  share(t,r,"l",s_reg,"extant")$denom(t,r,s_reg) = (1+tl(t,r))*pl.l(t,r)
                                                   *ld_ex.l(t,r,s_reg)*ld_ex0(r,s_reg)
                                                   /denom(t,r,s_reg);
  share(t,r,s,s_reg,"extant")$denom(t,r,s_reg) = pa.l(t,r,s)
                                                 *id_ex.l(t,r,s,s_reg)*id_ex0(r,s,s_reg)
                                                 /denom(t,r,s_reg);

else

* rebalance the cost shares in the case where any of the inputs
* with non-zero shares doesn't have associated inputs in the sector
  denom(t,r,s_reg) = sum(s, prod_share(s)$id0(r,s,s_reg))
                    +prod_share("k")$kd0(r,s_reg)
                    +prod_share("l")$ld0(r,s_reg);

* share of the shock for production with new capital
  share(t,r,s,s_reg,"new")$denom(t,r,s_reg) = prod_share(s)/denom(t,r,s_reg);
  share(t,r,"l",s_reg,"new")$denom(t,r,s_reg) = prod_share("l")/denom(t,r,s_reg);
  share(t,r,"k",s_reg,"new")$denom(t,r,s_reg) = prod_share("k")/denom(t,r,s_reg);

* share of the shock for production with extant capital
  share(t,r,s,s_reg,"extant")$denom(t,r,s_reg) = prod_share(s)/denom(t,r,s_reg);
  share(t,r,"l",s_reg,"extant")$denom(t,r,s_reg) = prod_share("l")/denom(t,r,s_reg);
  share(t,r,"k",s_reg,"extant")$denom(t,r,s_reg) = prod_share("k")/denom(t,r,s_reg);


);

* productivity shock for production with new capital
prod_ind(t,r,"k",s_reg,"new")$(kd.l(t,r,s_reg)*kd0(r,s_reg))
  = prod_ind(t,r,"k",s_reg,"new")*(1+share(t,r,"k",s_reg,"new")
                                     *shock(t,r,s_reg,"new")
                                     /((1+tk(t,r))*pr.l(t,r)*kd.l(t,r,s_reg)*kd0(r,s_reg)));
prod_ind(t,r,"l",s_reg,"new")$(ld.l(t,r,s_reg)*ld0(r,s_reg))
  = prod_ind(t,r,"l",s_reg,"new")*(1+share(t,r,"l",s_reg,"new")
                                     *shock(t,r,s_reg,"new")
                                     /((1+tl(t,r))*pl.l(t,r)*ld.l(t,r,s_reg)*ld0(r,s_reg)));
prod_ind(t,r,s,s_reg,"new")$(id.l(t,r,s,s_reg)*id0(r,s,s_reg))
  = prod_ind(t,r,s,s_reg,"new")*(1+share(t,r,s,s_reg,"new")
                                   *shock(t,r,s_reg,"new")
                                   /(pa.l(t,r,s)*id.l(t,r,s,s_reg)*id0(r,s,s_reg)));

* productivity shock for production with extant capital
prod_ind(t,r,"k",s_reg,"extant")$(kd_ex.l(t,r,s_reg)*kd_ex0(r,s_reg))
  = prod_ind(t,r,"k",s_reg,"extant")*(1+share(t,r,"k",s_reg,"extant")
                                        *shock(t,r,s_reg,"extant")
                                        /((1+tk(t,r))*pr_ex.l(t,r,s_reg)*kd_ex.l(t,r,s_reg)*kd_ex0(r,s_reg)));
prod_ind(t,r,"l",s_reg,"extant")$(kd_ex0(r,s_reg) and
                                    ld_ex.l(t,r,s_reg)*ld_ex0(r,s_reg))
  = prod_ind(t,r,"l",s_reg,"extant")*(1+share(t,r,"l",s_reg,"extant")
                                        *shock(t,r,s_reg,"extant")
                                        /((1+tl(t,r))*pl.l(t,r)*ld_ex.l(t,r,s_reg)*ld_ex0(r,s_reg)));
prod_ind(t,r,s,s_reg,"extant")$(kd_ex0(r,s_reg) and
                                  id_ex.l(t,r,s,s_reg)*id_ex0(r,s,s_reg))
  = prod_ind(t,r,s,s_reg,"extant")*(1+share(t,r,s,s_reg,"extant")
                                      *shock(t,r,s_reg,"extant")
                                      /(pa.l(t,r,s)*id_ex.l(t,r,s,s_reg)*id_ex0(r,s,s_reg)));
